!     Copyright (C) 2010 Benjamin Piaud
!
!     LIMBES is free software; you can redistribute it and/or modify
!     it under the terms of the GNU General Public License as published by
!     the Free Software Foundation; either version 3, or (at your option)
!     any later version.

!     LIMBES is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!     GNU General Public License for more details.

!     You should have received a copy of the GNU General Public License
!     along with LIMBES; if not, see <http://www.gnu.org/licenses/>



!      subroutine LIMBES_collision computes at the node (xi,yi)
!      the collision term (i.e. bgk for the moment):
!      dfdt = - (f-feq)/tau
!      first order euler scheme is used for temporal integration
!      LIMBES_f(1:LIMBES_nmax,1:LIMBES_nmax, 1:LIMBES_nvmax,1:LIMBES_nvmax)
!      -> f (real*8)
!      LIMBES_fnew(1:LIMBES_nmax,1:LIMBES_nmax, 1:LIMBES_nvmax,1:LIMBES_nvmax)
!      -> updated f (real*8)
!      LIMBES_vgradf(1:LIMBES_nmax,1:LIMBES_nmax, 1:LIMBES_nvmax,1:LIMBES_nvmax)
!      -> v.grad(f) (real*8)
!      needed argument
!      integer*8 xi,yi node coordinate

subroutine LIMBES_collision(ni)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: vi,ni
  real*8 :: tau,LIMBES_tau
  real*8 :: rho,ux,uy,RT,vx,vy,ax,ay
 
  !      computation of tau
  tau=LIMBES_tau(ni)

  !      computation of feq
  rho=LIMBES_rho(ni)
  ux=LIMBES_ux(ni)
  uy=LIMBES_uy(ni)
  RT=LIMBES_RT(ni)  

  !      forcing term
  call LIMBES_force(ni,ax,ay)

  ! collision step
  LIMBES_coll(ni,:)=-(LIMBES_f(ni,:)-LIMBES_feq(ni,:))/tau + (ax*(LIMBES_vx(:)-ux)/RT + ay*(LIMBES_vy(:)-uy)/RT)*LIMBES_feq(ni,:)

end subroutine LIMBES_collision

!----------------------------------------------------------------
!      function LIMBES_tau computes the collision time tau
!      which can depends on rho and RT
!      needed argument
!      integer*8 xi,yi node coordinate

real*8 function LIMBES_tau(ni)
  use LIMBES_mod_var
  use LIMBES_mod_fluid_config
  implicit none
  integer*8 :: ni
  real*8 :: tau

  select case (FLUID_tau_expr)
  case default
     print*,'error from LIMBES_tau: tau exression not known'
     print*,'FLUID_tau_expr=',FLUID_tau_expr
     stop

  case ('constant_tau')
     tau=FLUID_tau

  end select

  LIMBES_tau=tau
end function LIMBES_tau
